Winding up by a quench: vortices in the wake of rapid Bose-Einstein condensation 
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A second order phase transition induced by a rapid quench 
can lock out topological defects with densities far exceeding 
their equilibrium expectation values. We use quantum kinetic 
theory to show that this mechanism, originally postulated in 
the cosmological context, and analysed so far only on the 
mean field classical level, should allow spontaneous generation 
of vortex lines in trapped Bose-Einstein condensates of simple 
topology, or of winding number in toroidal condensates. 
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An as yet unachieved goal of experiments on trapped 
ultra-cold alkali gases |TJ is the exhibition of a persistent 
vortex. Since the reason that superfluid vortices are per- 
sistent is that there is a high energetic barrier between 
the metastable vortex state and the non-rotating true 
ground state, spinning up a non-rotating condensate once 
it is fully grown seems likely to be difficult to accomplish 
without excessive heating. In this Letter we show that 
a rotating condensate may instead grow spontaneously 
from fluctuations during a non-equilibrium quench. Not 
only is this a possible procedure for generating vortices: it 
actually provides a strong additional motivation for vor- 
tex experiments, by making vortex production in cooled 
alkali gases a test of a fundamental prediction of non- 
equilibrium statistical mechanics. 

The widely applied time-dependent Ginzburg-Landau 
theory (TDGL) predicts failure of equilibrium during a 
second order phase transition at finite speed. In TDGL, 
the complex order parameter ip(r,t) obeys 



W = /3(^-V 2 + M -Ah/f)<A, 



(1) 



where = (fcflT) -1 , and To and A > are phenomeno- 
logical parameters. The thermodynamical variable /i be- 
haves near the critical point, in the case we consider, as 



/i= -(T c -T)+0(T c -T) 2 , 



(2) 



where T c is the critical temperature. The equilibration 
time for long wavelengths is r = ToksT '/\fi\. The sys- 
tem's disordered phase is described by (i < 0, so that 
ijj = is a stable fixed point of (m. The ordered phase 
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appears when \x > 0, since then the stable fixed points lie 
on the circle \ip\ 2 = /it/A, and the phase 9 of tp = \ifj\e 
becomes a new macroscopic variable. 

A quench occurs if [i changes with time from negative 
to positive values. The divergence of the equilibration 
time r at the critical point [i = is associated with crit- 
ical slowing down. Because of this critical slowing down, 
^ //i must exceed I /r in some neighbourhood of the crit- 
ical point, and so there must be an epoch in which the 
system is out of equilibrium. What are at the beginning 
of this epoch mere fluctuations in the disordered phase, 
in which higher energy modes happen momentarily to 
be more populated than the lowest mode, can thus pass 
unsuppressed by equilibration into the ordered phase, to 
become topologically non-trivial configurations of the or- 
der parameter field ip. One therefore expects topological 
defects, such as vortex lines, to form spontaneously dur- 
ing a transition at sufficient speed [|| . 

The interval within which equilibration is negligible 
can be identified as the period wherein \t\/r < 1. If we 
define the quench time scale tq by letting = </tq 
(choosing t — as the moment the system crosses the 
critical point), this implies that the crucial interval is 
—i<t<t, for i = y/TQTo The correlation length 

£ for fluctuations at time t = — t is then given by 
h/(2M£ 2 ) = which (assuming T(-f) = T c ) im- 

plies that | = A Tc (tq/t- ) 1/4 , for A T = h(2Mk B T)-^ the 
thermal de Broglie wavelength. Taking this correlation 
length as giving the typical domain size surrounding a 
defect H implies that the vortex line density, in bulk, 
should be proportional to Tq 1 ^ 2 [Q. Alternatively one 
can consider the transition to occur within a toroidal ves- 
sel, so that independent random settings of the order pa- 
rameter phase, at different points around the torus, can 
produce a net vorticity, W = =- §dl V#. This implies a 
superflow, with velocity UV6/M Q. In this case one es- 
timates one independently chosen phase within each cor- 
relation length £; modeling the phase distribution around 
the torus as a random walk suggests that the net vorticity 
should be proportional to £~ 1 ^ 2 , hence to Tq 1 ^ 8 j2|. 

Although ingenious experiments have recently been 
performed to test this theory, in liquid helium JE|], and 
numerical studies have supported its scaling predictions 
Q , it would be even more informative to have analogous 
results in a weakly interacting system, such as a dilute 
trapped alkali gas. Assuming tq is the scattering time, 
evaporative cooling techniques yield (tq/tq) 1 ^ 4 of order 
one, and so £ is essentially At c ■ For atoms at several hun- 
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dred nK, this means £ ~ 100 nm, smaller than current 
condensates. As numerical simulations show ||, this is a 
generously low lower bound on the distance between vor- 
tex lines, but it does indicate that spontaneous vorticity 
should be within experimental reach. Considering this 
intriguing prospect raises an obvious question: is TDGL 
actually relevant to finite samples of dilute gas, far from 
equilibrium? 

We therefore begin again from first principles, and con- 
sider a dilute Bose gas in a trap, with the Hamiltonian 



H 



(5) 



where n k (t) = Tt(ph k ). The equation governing the 
n k (t) follows simply from (Eh: 



h k = r e^[l + (1 - e^-^fifc] 



(6) 



This equation may be integrated for general /3(t),p(t). 
But as a simple form valid near the critical point, we im- 
pose (3(t)[fj,(t) — E k ] = {t — &k)/TQ, defining tq as well as 
the bias time scales The n k {t) that result, from the 



— / d 3 r (JV?/'! 2 + U {r)^ijj + ina^ 2 ^ 2 ) (3) equilibrium initial values n k {U) 



•1 



where ip(r) annihilates a boson at position r, U gives the 
trap potential, and a is the s-wave scattering length. As 
always, ip{r) = 'Yl k U k {r)'ip k defines a decomposition of 
the system into orthogonal modes described by single- 
particle wave functions u k - In the earliest stages of con- 
densation, it is sufficient to take the single-particle energy 
eigenstates as defining the normal modes of the gas. 

We now construct a quantum kinetic theory (QKT), 
by considering the lowest energy modes of the trap, up 
to some energy En, to be an open quantum system (the 
'condensate band'), interacting via two-particle s-wave 
scattering with the higher modes, treated as a 'reservoir 
band' Q. We model evaporative cooling by prescrib- 
ing that the reservoir band is always in equilibrium, but 
with a time dependent temperature f3~ l {t) and chem- 
ical potential p(t) (which can become positive as long 
as it remains below Er). We then form the reduced 
density operator for the condensate band by tracing out 
the reservoir. The condensate band will not remain in 
equilibrium with the reservoir; the time evolution of its 
reduced density operator is the problem to be solved. 

In the earliest stages of condensation, before nonlinear 
coherent interactions become important, one can derive a 
simple master equation for the condensate band, strongly 
reminiscent of that of a multi-mode laser: 



' \ ih 



(h k p + pn k ) - p J 



(4) 



where E k are the energies of the normal modes. The Tk 
are scattering rates, which may be computed; they will 
generally be of the order of the Boltzmann scattering 
rate. We actually expect the /c-dependence of the Tk to 
be weak as long as the temperature is much larger than 
the trap level spacing, so we will hereafter replace I\. with 
Tq, which will play exactly the same role as 1/tq did in 
TDGL. The non-Hermitian part of (EJ) is due to collisions 
in which one particle leaves or joins the condensate for 
or from the reservoir. 

An ansatz which solves (0) is furnished by 



are incomplete Gamma functions; they only begin to de- 
part significantly from their instantaneous equilibrium 
values after t — ~ — iJtq/To = —t. Past these points, 
the fik lag below their equilibrium values. This clarifies 
the effect of the critical slowing down: as Bose enhance- 
ment turns on, the rates of scattering into the condensate 
increase; but the numbers of particles required by equilib- 
rium increase faster still, and so the ability of scattering 
to maintain equilibrium rapidly declines. 



After these times, we can approximate (1 — e 



0(E k -ii) 



(t — /tq and match to equilibrium at early times, to 
see that 



n k {t) =T e2 



dt e 



(7) 



For times after t—'&k —i> each fik grows explosively, be- 
cause the atomic scattering analogue of stimulated emis- 
sion into the fcth mode is turning on strongly: n k is be- 
coming large enough that the term proportional to it on 
the RHS of (||) dominates the other term. Bose-enhanced 
scattering then enables the mode to begin a very rapid 
'whiplash' to catch up with equilibrium. So the interval 
— t < t < $k + t is indeed a transition zone between 
equilibrium above T c , and the onset of coherent processes 
below T c . It is obvious that for a higher energy mode to 
have any significant chance of competing successfully for 
particles with the lowest mode, it cannot afford to begin 
explosive growth much later than the lowest mode. This 
implies that "Dk < t, or /3E k < (TqTq)" 1 / 2 , limits the 
range of significantly competitive modes. Since in bulk 
or in a toroidal trap we have E k oc k 2 , this gives 



i=l = h(2Mk B T c )-?{T Q T Q ) 
k 



1/4 



(8) 



which is the same conclusion reached by TDGL g. 

For the toroidal problem, the density operator pre- 
scribed by the linear quantum kinetic theory is equiv- 
alent to a distribution of coherent states with probabili- 
ties proportional to exp — ^ fe ^\i^ k \ 2 , for Fourier modes 
k. While W is not a simple function of xj) k , the idea 
that there are as many independent random phases as 
non-negligible fik (i) still seems reasonable, and we expect 
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typical vorticities of order y fcL/27r, for L the perimeter 
of the torus. This again co-incides with the TDGL pre- 
diction. 

Our conclusion at this point is that QKT agrees with 
the phenomenological theory, in predicting that for suf- 
ficiently rapid quenches the probability of forming a 
small 'seed' of condensate with non-zero vorticity is of 
order one. But since superfluid currents only become 
metastable above a threshold condensate density, not 
all of this initial vorticity will survive as the conden- 
sate grows. To follow the non-equilibrium evolution of 
a trapped condensate into the non-linear regime, with 
quantum kinetic theory, is a challenging problem. We 
therefore restrict our analysis to a simple toy model, 
which affords some qualitative insight, and allows a com- 
parison between TDGL and QKT. 

The toy model replaces the condensate band of many 
low energy modes by a system with only two modes, 
representing states with two different angular momenta. 
Because the self-Hamiltonian for this two-mode system 
must conserve both particle number and angular momen- 
tum, it must conserve separately the numbers of particles 
in both modes. We therefore choose 
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H = E[nx + yy(ni + h{ + 4hih )] 



(9) 



Because we have incorporated the Bose enhancement of 
inter-mode repulsion (the factor of 4 instead of 2 in front 
of the h\hi term, which is of course the best case value, 
obtained when uq and u\ overlap completely), we make 
the state with all particles in the 1 mode a local minimum 
of the energy for n 1+77,2 > (N c +1). For two lowest modes 
of a typical oblate magneto-optical trap, we have (3E of 
order 10 -2 ; for proposed toroidal traps with perimeter of 
order 10 -2 cm, at similar temperatures, (3E could be as 
low as 10~ 5 . (Rotating the gas before condensation could 
also lower the effective energy bias, and even favour ro- 
tating states over the ground state.) The experimental 
range of N c is around 100 for compact traps, but as low 
as 1 for the torus; this does not take into account the 
Thomas-Fermi expansion of the condensate wave func- 
tion, which in fact can make N c rise significantly at large 
particle numbers. 

We also assume interactions between both condensate 
modes and the quasi-continuum of reservoir modes, of 
the form implied by the Hamiltonian (|^). Upon tracing 
over the dilute gas reservoir, we obtain a master equa- 
tion of more complicated form than (|j), which includes 
saturation effects, as well as scattering of reservoir atoms 
off the condensate (with no resulting change in the con- 
densate number). For present purposes only the diagonal 
part of this equation is necessary: 

PriQ,ni — r(^) [Rno,ni Rno — l.ni ~t~ 5Vio,ni ^Tlo,Tli — l] 

— r(t)[T no+ i jni — T n0: , ll+1 ] 



R n0ini = (n Q + l)[e 0fi Pnoini - e- {no+2ni) Pno+hni ] 
S no , ni = (m + l)[e^p„ , ni - e^ ( ^ +2 " 0+ni) p„ ,„ 1+1 ] 

Pri — l,ni 

Pn , ni -l\ , (10) 



where T(t) and f (t) are again scattering rates (for scat- 
tering into/out of the condensate, and off the condensate, 
respectively) which may be computed for any specific 
condensate-reservoir coupling. Wc will hereafter assume 
r = (3ET, which is accurate for simple trap configura- 
tions when the temperature is much larger than the trap 
level spacing. (This f3E factor justified neglecting these 
bouncing-off processes in the linear regime; it appears be- 
cause most reservoir particles are so much faster than the 
condensate particles that they are unlikely to strike them 
without dislodging them from the condensate band.) 

Equation ( |io| ) provides a complete description of con- 
densation in the toy model, including initial seeding from 
fluctuations, coherent growth, relaxation into metastable 
states, and eventual equilibration by thermal barrier 
crossing. While it would be straightforward to solve 
numerically, we can obtain more understanding of the 
growth process by extracting from it an equation of mo- 
tion for no and n\ . This may be done, among other ways, 
by taking uq — > Nx and n\ — > Ny for continuous x and y 
and N of order (/3i?) _1 . Expanding the finite differences 
in (|l0|) in powers of derivatives with respect to x and y, 
one obtains a Fokker-Planck-likc equation, the Liouville 
terms of which describe a flow along deterministic tra- 
jectories in (x, y)-space. Dropping higher order terms in 
1/iV (since these are significant only at small no, n>i, when 
diffusion dominates systematic evolution but we are able 
to use the linear analysis described above), these trajec- 
tories obey 



h = Tn 



_ ^(«0+2ni) 



+2flEme~2 



m.\N c +no-m\ • u P E 



sinh (N c + no - ni) 



h± = Tn-i 



(11) 



-2f3Enoe~-^ N c +n °- ni1 S mh^-(N C + n - m) 



The first question is, how important is the systematic 



evolution prescribed by (11) compared to the diffusive 
evolution also contained in (|l(])? We can address this 
question by examining a Gaussian approximation to the 
Fokker-Planck equation from which ( |TT| ) came. Fig. 1 
shows selected solutions to ([n]) together with 68% proba- 
bility contours for Gaussian approximations to p(no, ni), 
starting from initial delta functions. Fig. 1(a) shows that 
for a slow quench, diffusion is in fact very strong; this 
does not necessarily mean that the metastable state is 
not reached, but that it may be reached by diffusive nu- 
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cleation rather than via the critical slowing down mech- 
anism we are considering. Even here, though, there are 
'channels' near the axes in which diffusion is weaker. For 
fast quenches, as shown in Fig. 1(b), diffusion is clearly a 
small correction to predominantly systematic evolution. 
In such cases, therefore, we may obtain accurate esti- 
mates of the probability of reaching the metastable state, 
by using the linear analysis described above to compute 
the distribution p(no,m) at some 'coherent start time' 
t s ~ t, and then letting the distribution flow under (11). 
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FIG. 1. Some solutions to (|ll|), with diffusion illustrated by 
68% probability contours at selected times. Axes are numbers 
of particles in modes (horizontal) and 1 (vertical) ; starting 
points are on line no + ni = (2f3E)~ 1 , with no(t s ) = (2/3E)~ 1 
setting the start time t s . Quench is j3fj, — tanh(t/rQ). 
N c = 10; other parameters are (a) f3E = 10~ 2 , Ttq = 40, 
and (b) (3E = 10~ 3 , Ttq = 10. 

Having established that the systematic evolution of 
( |lT| ) provides a good description, after i, of a fast quench 
in the toy model, we can now compare it to the TDGL 
evolution. When uq + 2ni and N c + ni + 2uq are both 
close to N c n/E, or for low enough particle numbers, 
the first line in each equation of ( pi] ) is indeed equiva- 
lent to a TDGL equation (as may be seen by replacing 
rij — > IV'jl 2 )- But the second line in each equation is 
not of Ginzburg-Landau form: it does not involve «, and 
the expression it implies for tpj is not a gradient with re- 
spect to tp* . These non-GL terms conserve uq + ni, and 
describe doubly Bose-enhanced dissipation due to scat- 
tering of reservoir particles off the condensate. They turn 
out to imply that the system equilibrates in energy faster 
than it equilibrates in particle number. 

Some representative solutions to ([ll]) are shown in 
Fig. 2, together with the |?/>j| 2 given by the TDGL 
equation. It is clear that for sufficiently fast quenches, 
the two theories accord quite well, but that for slower 
quenches TDGL significantly overestimates the proba- 
bility of reaching the metastable state. If our two modes 
are taken to be different Fourier modes in a toroidal trap, 
the vorticity of a state is simply the vorticity of the more 
populated mode, so that the line no — n\ is the border 
between vorticities; all initial points in Fig. 2 are above 
this line. So not even TDGL evolution conserves vortic- 
ity, but the QKT evolution changes vorticity more easily, 
especially for slower quenches. 
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FIG. 2. Trajectories from QKT (solid) and TDGL (dot- 
ted); heavy dashed line is threshold for metastability of 
mode 1. Initial times are t; quench is /3fi — tanh(t/rQ), 
P = /3 c e tanh(t/T o ) . Parameters are N c = 100, and (a) 
Ttq = 10, f3 c E = 0.01; (b) Ttq = 100, (5 C E = 0.05. 

Despite the shortcomings of TDGL revealed by our 
toy model, we would like to emphasize that in fact QKT 
does show that TDGL is relevant to trapped dilute gases, 
even very far from equilibrium: what TDGL requires is 
not outright rejection, but corrections, from diffusion and 
dissipation. And although these corrections may be sub- 
stantial, the gross features predicted by TDGL are still 
recovered, with faster quenches and smaller biases. While 
the extension of quantum kinetic theory beyond toy mod- 
els, to realistic descriptions of topological defect forma- 
tion, will obviously require much further study, we believe 
that the prospects for experimental realization of spon- 
taneous defects, as predicted by the Ginzburg-Landau 
theory, are very encouraging. 
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